#' ---
#' title: "Replication script for Reporting the polls - The quality of media reporting of vote intention polls in the Netherlands"
#' author: "Tom Louwerse and Rozemarijn van Dijk"
#' ---


# Load packages -----------------------------------------------------------
library(tidyverse)
library(lubridate)
library(rio)
library(effects)
library(lme4)
library(texreg)
library(date)
library(ggplot2)
library(estimatr)
library(miceadds)
library(sandwich)


# Load data ---------------------------------------------------------------

# Data on polls
data_polls <- readRDS("data_polls.rds") %>%
  mutate(year = format(Article_date, "%Y")) %>%
  filter(year %in% c(2010, 2012, 2017)) %>%
  mutate(Poll_publication_date = as_datetime(Poll_publication_date))

# Information on survey report quality
data_polls_moe <- readRDS("data_pollsters.rds") 

# Information on media
data_media <- readRDS("data_media_type.rds")

# Merge data
data <- as_tibble(data_polls) %>%
  left_join(data_polls_moe, 
                  by=c("Poll_company"="Poll_company", 
                       "Poll_publication_date"="Poll_publication_date")) %>%
  left_join(data_media, 
                  by=c("Article_source")) %>%
  mutate(TV=case_when(Tabloid == 1 ~ 0,
                      Broadsheet == 1 ~ 0,
                      TRUE ~ 1))


# Descriptive data on WAPOR items -----------------------------------------


# WAPOR score per year per article
count_WAPOR <- data %>%
  filter(Article_poll_type %in% c(2, 3)) %>%
  group_by(Article_ID) %>%
  summarize(WAPOR=max(WAPOR, na.rm = TRUE), year=mean(as.numeric(year), na.rm=TRUE)) %>%
  ungroup() %>%
  group_by(year, WAPOR) %>%
  summarize(N_articles=length(Article_ID))


# WAPOR score per year per article
count_WAPOR_items <- data %>%
  filter(!is.na(Poll_company)) %>%
  filter(Article_poll_type %in% c(3)) %>%
  group_by(Article_ID) %>%
  summarize_at(vars(Poll_universe, `Poll_sampling method`,Poll_non_response_rate,Poll_sample_size,
                    Poll_weights, Poll_dates, Poll_data_collection_method,
                    Margin_of_error, year),
               max, na.rm=TRUE) %>%
  ungroup() %>%
  group_by(year) %>%
  summarise(across(!Article_ID, .fns=function(x) round(mean(x) * 100, 1)))

rownames(count_WAPOR_items) <- count_WAPOR_items$year 
count_WAPOR_items <- count_WAPOR_items %>%
  select(-year) %>%
  t() %>%
  xtable::xtable(digits=1)

count_WAPOR_items
print(count_WAPOR_items, type="html", file="desc_table_articles.doc", 
      html.table.attributes="border=0")

count_WAPOR_items_n_articles <- data %>%
  filter(!is.na(Poll_company)) %>%
  filter(Article_poll_type %in% c(3)) %>%
  group_by(Article_ID) %>%
  summarize_at(vars(Margin_of_error, year),
               max, na.rm=TRUE) %>%
  ungroup() %>%
  group_by(year) %>%
  summarise(across(!Article_ID, .fns=list(Number_articles=length)))
count_WAPOR_items_n_articles


count_WAPOR_items_exclPW <- data %>%
  filter(Poll_company!="Peilingwijzer") %>%
  filter(Article_poll_type %in% c(3)) %>%
  group_by(Article_ID) %>%
  summarize_at(vars(Poll_universe, `Poll_sampling method`,Poll_non_response_rate,Poll_sample_size,
                    Poll_weights, Poll_dates, Poll_data_collection_method,
                    Margin_of_error, year),
               max, na.rm=TRUE) %>%
  ungroup() %>%
  group_by(year) %>%
  summarise(across(!Article_ID, .fns=function(x) round(mean(x) * 100, 1)))
count_WAPOR_items_exclPW


rownames(count_WAPOR_items_exclPW) <- count_WAPOR_items_exclPW$year 
count_WAPOR_items_exclPW <- count_WAPOR_items_exclPW %>%
  select(-year) %>%
  t() %>%
  xtable::xtable(digits=1)


print(count_WAPOR_items_exclPW, type="html", file="desc_table_articles_exclpw.doc", 
      html.table.attributes="border=0")


count_WAPOR_items_pollsters <- data_polls_moe %>%
  mutate(year = format(Poll_publication_date, "%Y")) %>%
  group_by(year) %>%
  select(Pollster_poll_universe, `Pollster_poll_sampling method`, Pollster_poll_non_response_rate,
         Pollster_poll_sample_size, Pollster_poll_weights, Pollster_poll_dates, Pollster_poll_data_collection_method,
         Pollster_reported_MOE) %>%
  summarize_all(mean) %>%
  mutate_at(vars(-year), ~ round(. * 100, 1))

rownames(count_WAPOR_items_pollsters) <- count_WAPOR_items_pollsters$year 
count_WAPOR_items_pollsters <- count_WAPOR_items_pollsters %>%
  select(-year) %>%
  t() %>%
  xtable::xtable(digits=1)

count_WAPOR_items_pollsters
print(count_WAPOR_items_pollsters, type="html", file="desc_table_pollsters.doc",
      html.table.attributes="border=0")


# Number of unique articles in analysis (cat 2 or 3)
data %>% 
  filter(Poll_company %in% c("LISS", "Peil", "Kantar", "I&O Research", "Peilingwijzer", "Ipsos/1V", "EenVandaag")) %>%
  filter(!is.na(WAPOR), !is.na(Margin_of_error), !is.na(Poll_publication_date), !is.na(Poll_company), Poll_company != "NA") %>%
  filter(Article_poll_type %in% c(2,3)) %>%
  filter(Party_result_explicit ==1 | Article_poll_type == 3) %>%
  select(Article_ID) %>%
  distinct() %>%
  nrow()
# 392

# Number of unique polling-companies in articles
data %>% 
  filter(Poll_company %in% c("LISS", "Peil", "Kantar", "I&O Research", "Peilingwijzer", "Ipsos/1V", "EenVandaag")) %>%
  filter(!is.na(WAPOR), !is.na(Margin_of_error), !is.na(Poll_publication_date), !is.na(Poll_company), Poll_company != "NA") %>%
  filter(Article_poll_type %in% c(2,3)) %>%
  filter(Party_result_explicit ==1 | Article_poll_type == 3) %>%
  select(Article_ID, Poll_company) %>%
  distinct() %>%
  nrow()
# 524


# H1: The higher the number of WAPOR items included in a pollster  --------

# Create data frame for this analysis
# This aggregates the data to the article-pollster level (i.e. a mention of a polling company in an article)

data_h1 <- data %>%
  filter(Poll_company %in% c("LISS", "Peil", "Kantar", "I&O Research", "Peilingwijzer", "Ipsos/1V", "EenVandaag")) %>%
  filter(!is.na(WAPOR), !is.na(Margin_of_error), !is.na(Poll_publication_date), !is.na(Poll_company), Poll_company != "NA") %>%
  filter(Article_poll_type %in% c(2,3)) %>%
  filter(Party_result_explicit ==1 | Article_poll_type == 3) %>%
  mutate(PW = as.numeric(Poll_company == "Peilingwijzer")) %>%
  mutate(NOS = as.numeric(as.numeric(Article_source == "NOS"))) %>%
  mutate(year = as.factor(year)) %>%
  group_by(Article_ID, Poll_company, year) %>%
  summarise(WAPOR=mean(WAPOR),
            Poll_universe = max(Poll_universe, na.rm=TRUE),
            Poll_sampling_method = max(`Poll_sampling method`, na.rm=TRUE),
            Poll_non_response_rate = max(Poll_non_response_rate, na.rm=TRUE),
            Poll_sample_size = max(Poll_sample_size, na.rm=TRUE),
            Poll_dates = max(Poll_dates, na.rm=TRUE),
            Poll_data_collection_method = max(Poll_data_collection_method, na.rm=TRUE),
            Margin_of_error=as.numeric(mean(Margin_of_error)>0.5),
            PW=as.factor(as.numeric(mean(PW, na.rm=TRUE)>0.5)),
            Tabloid = max(Tabloid, na.rm=TRUE),
            Broadsheet = max(Broadsheet, na.rm=TRUE),
            TV = max(TV, na.rm=TRUE),
            NOS=as.factor(as.numeric(mean(NOS, na.rm=TRUE)>0.5)),
            Article_source = as.factor(first(Article_source)), 
            Article_length = mean(Article_length, na.rm=TRUE),
            Poll_company = factor(first(Poll_company)),
            Article_poll_type = factor(first(Article_poll_type), levels=c(2,3)),
            Pollster_WAPOR_index = mean(Pollster_WAPOR_index),
            Pollster_reported_MOE = factor(as.numeric(mean(Pollster_reported_MOE, na.rm=TRUE)>0.5), 
                                           levels=c(0,1)),
            Pollster_poll_universe = factor(as.numeric(mean(Pollster_poll_universe, na.rm=TRUE)>0.5), 
                                            levels=c(0,1)),
            Pollster_poll_sampling_method = factor(as.numeric(mean(`Pollster_poll_sampling method`, na.rm=TRUE)>0.5), 
                                                   levels=c(0,1)),
            Pollster_poll_non_response_rate = factor(as.numeric(mean(Pollster_poll_non_response_rate, na.rm=TRUE)>0.5), 
                                                     levels=c(0,1)),
            Pollster_poll_sample_size = factor(as.numeric(mean(Pollster_poll_sample_size, na.rm=TRUE)>0.5), 
                                               levels=c(0,1)),
            Pollster_poll_dates = factor(as.numeric(mean(Pollster_poll_dates, na.rm=TRUE)>0.5), 
                                         levels=c(0,1)),
            Pollster_poll_data_collection_method = factor(as.numeric(mean(Pollster_poll_data_collection_method, na.rm=TRUE)>0.5), 
                                                          levels=c(0,1))) %>%
  ungroup() %>%
  mutate(year=factor(year)) %>%
  mutate(type_of_moe = factor(case_when(Pollster_reported_MOE=="0" ~ "No MOE",
                                        Pollster_reported_MOE=="1" & Poll_company !="Peilingwijzer" ~ "MOE",
                                        Pollster_reported_MOE=="1" & Poll_company =="Peilingwijzer" ~ "Range"),
                              levels=c("No MOE", "MOE", "Range")))



# Model H1: WAPOR items
model100 <- lmer(WAPOR ~ Pollster_WAPOR_index + log(Article_length) + Article_poll_type + year + (1| Poll_company), 
                 data=data_h1 %>% filter(Poll_company != "Peilingwijzer"))
summary(model100)
plot(allEffects(model100))

# Clustered standard errors with year fixed effects
model100c <- lm.cluster(WAPOR ~ Pollster_WAPOR_index + log(Article_length) + Article_poll_type + year, 
                        data=data_h1, 
                        cluster='Poll_company',
                        subset=data_h1$Poll_company!="Peilingwijzer")
summary(model100c)


## WAPOR item models - separate

# Sample size (appendix only)
model101 <- glmer(Poll_sample_size ~ Pollster_poll_sample_size + log(Article_length) + Article_poll_type + year + (1| Poll_company), data=data_h1,
                  subset=data_h1$Poll_company!="Peilingwijzer", 
                  family=binomial(link="logit"))
summary(model101)


model101c <- glm.cluster(Poll_sample_size ~ Pollster_poll_sample_size + log(Article_length) + Article_poll_type + year, 
                         data=data_h1, 
                         cluster='Poll_company',
                         subset=data_h1$Poll_company!="Peilingwijzer", 
                         family=binomial(link="logit"))
summary(model101c)



# Data collection method
table(data_h1$Poll_data_collection_method, data_h1$Pollster_poll_data_collection_method)

model102 <- glmer(Poll_data_collection_method ~ Pollster_poll_data_collection_method + log(Article_length) + Article_poll_type + year + (1| Poll_company), data=data_h1,
                  subset=data_h1$Poll_company!="Peilingwijzer", 
                  family=binomial(link="logit"))
summary(model102)
plot(allEffects(model102))


model102c <- glm.cluster(Poll_data_collection_method ~ Pollster_poll_data_collection_method + log(Article_length) + Article_poll_type + year, 
                         data=data_h1, 
                         cluster='Poll_company',
                         subset=data_h1$Poll_company!="Peilingwijzer", 
                         family=binomial(link="logit"))
summary(model102c)



# H2: Pollster MOE  --> News report MOE -----------------------------------

model200 <- glmer(Margin_of_error ~ Pollster_reported_MOE + log(Article_length) + Article_poll_type +
                    year + (1| Poll_company), data=data_h1, family=binomial(link="logit"))
summary(model200)
plot(allEffects(model200))


# Fixed effects version
model200c <- glm.cluster(Margin_of_error ~ Pollster_reported_MOE + log(Article_length) + Article_poll_type + year, 
                         data=data_h1, 
                         cluster='Poll_company',
                         family=binomial(link="logit"))

summary(model200c)


# H3: Pollster MOE type  --> News report MOE ----------------------------------
model300 <- glmer(Margin_of_error ~ type_of_moe + log(Article_length) + Article_poll_type +
                    year + (1| Poll_company), 
                  data=data_h1, 
                  family=binomial(link="logit"))
summary(model300)
plot(allEffects(model300))


# Fixed effects version
model300c <- glm.cluster(Margin_of_error ~  type_of_moe + log(Article_length) + Article_poll_type + year, 
                         data=data_h1, 
                         cluster='Poll_company',
                         family=binomial(link="logit"))

summary(model300c)



# H4A: Pollster MOE --> Difference ---------------------------------------------

data_h4a <- data %>%
  filter(!is.na(WAPOR), !is.na(Margin_of_error), !is.na(Poll_company), Poll_company != "NA") %>%
  filter(Article_poll_type %in% c(2, 3)) %>%
  mutate(PW = as.factor(as.numeric(Poll_company == "Peilingwijzer"))) %>%
  mutate(NOS = as.factor(as.numeric(Article_source == "NOS"))) %>%
  mutate(year = as.factor(year)) %>%
  mutate(Article_poll_type = factor(Article_poll_type)) %>%
  filter(!is.na(Report_party_difference_significant)) %>%
  mutate(Difference_correct = as.numeric(Report_party_difference_significant==1 | Report_party_difference_significance_discussed == 1)) %>%
  filter(!is.na(Difference_correct)) %>%
  mutate(year2017=as.factor(year==2017)) %>%
  mutate(type_of_moe = factor(case_when(Pollster_reported_MOE=="0" ~ "No MOE",
                                        Pollster_reported_MOE=="1" & Poll_company !="Peilingwijzer" ~ "MOE",
                                        Pollster_reported_MOE=="1" & Poll_company =="Peilingwijzer" ~ "Range"),
                              levels=c("No MOE", "MOE", "Range")))


model400 <- glmer(Difference_correct ~ type_of_moe + log(Article_length) + Article_poll_type + year + (1 | Article_ID), 
                  data=data_h4a, 
                  family=binomial(link="logit"),
                  control= glmerControl(optimizer="bobyqa", optCtrl = list(maxfun=1e7)))
summary(model400)
plot(allEffects(model400), type="response")

# Fixed effects specification
model400c <- glm.cluster(Difference_correct ~ type_of_moe + log(Article_length) + Article_poll_type + year, 
                         data=data_h4a, 
                         cluster='Article_ID',
                         family=binomial(link="logit"))
summary(model400c)


# Some additional descriptives
data_h4a %>%
  group_by(type_of_moe) %>%
  summarise(mean=mean(Difference_correct))

t.test(Difference_correct ~ PW, data=data_h4a)

table(data_h4a$Difference_correct) / sum(table(data_h4a$Difference_correct))



# H4B: Pollster MOE --> Change ---------------------------------------------
data_h4b <- data %>%
  filter(Report_change==1) %>%
  filter(!is.na(WAPOR), !is.na(Margin_of_error)) %>%
  filter(Article_poll_type %in% c(2, 3)) %>%
  filter(Poll_company != "NA") %>%
  mutate(PW = as.factor(as.numeric(Poll_company == "Peilingwijzer"))) %>%
  mutate(NOS = as.factor(as.numeric(Article_source == "NOS"))) %>%
  mutate(year = as.factor(year)) %>%
  mutate(Article_poll_type = factor(Article_poll_type, levels=c(2,3))) %>%
  filter(!is.na(Report_change_significant)) %>%
  mutate(Change_correct = as.numeric(Report_change_significant ==1 | Report_change_significance_discussed == 1)) %>%
  filter(!is.na(Change_correct)) %>%
  mutate(type_of_moe = factor(case_when(Pollster_reported_MOE=="0" ~ "No MOE",
                                        Pollster_reported_MOE=="1" & Poll_company !="Peilingwijzer" ~ "MOE",
                                        Pollster_reported_MOE=="1" & Poll_company =="Peilingwijzer" ~ "Range"),
                              levels=c("No MOE", "MOE", "Range")))





model451 <- glmer(Change_correct ~ type_of_moe + log(Article_length) + Article_poll_type + year +  (1 | Article_ID), 
                  data=data_h4b, family=binomial(link="logit"))
summary(model451)
plot(allEffects(model451))

# Fixed effects specification
model451c <- glm.cluster(Change_correct ~ type_of_moe + log(Article_length) + Article_poll_type + year, 
                         data=data_h4b, 
                         cluster='Article_ID',
                         family=binomial(link="logit"))
summary(model451c)


# Some additional descriptives
table(data_h4b$Change_correct) / sum(table(data_h4b$Change_correct)) # % of changes correctly reported


data_h4b %>%
  group_by(type_of_moe) %>%
  summarise(mean=mean(Change_correct))




# H5: Tabloid/broadsheet --------------------------------------------------
# These models replicate models for the previous hypotheses, but including variables Broadsheet & TV

# WAPOR
model501 <- lmer(WAPOR ~ Broadsheet + TV + Pollster_WAPOR_index + log(Article_length) + Article_poll_type + year + (1| Poll_company), data=data_h1,
                 subset=data_h1$Poll_company!="Peilingwijzer")
summary(model501)


# Fixed effects specification
model501c <- lm.cluster(WAPOR ~ Broadsheet + TV + Pollster_WAPOR_index + log(Article_length) + Article_poll_type + year, 
                        data=data_h1, 
                        cluster='Poll_company',
                        subset=data_h1$Poll_company!="Peilingwijzer")

summary(model501c)



# MOE
model502 <- glmer(Margin_of_error ~ Broadsheet + TV + type_of_moe + log(Article_length) + Article_poll_type +
                    year + (1| Poll_company), data=data_h1, family=binomial(link="logit"))
summary(model502)

# Fixed effects specification
model502c <- glm.cluster(Margin_of_error ~   Broadsheet + TV + type_of_moe + log(Article_length) + Article_poll_type + year, 
                         data=data_h1, 
                         cluster='Poll_company',
                         family=binomial(link="logit"))

summary(model502c)




# Difference
model503 <- glmer(Difference_correct ~ Broadsheet + TV + type_of_moe + log(Article_length) + Article_poll_type + year + (1 | Article_ID), 
                  data=data_h4a, 
                  family=binomial(link="logit"), 
                  control=glmerControl(optimizer="bobyqa", 
                                       optCtrl=list(maxfun=2e5)))
summary(model503)


# Fixed effects specification
model503c <- glm.cluster(Difference_correct ~ Broadsheet + TV + type_of_moe + + log(Article_length) + Article_poll_type + year, 
                         data=data_h4a, 
                         cluster='Article_ID',
                         family=binomial(link="logit"))
summary(model503c)



# Change
model504 <- glmer(Change_correct ~ Broadsheet + TV + type_of_moe + log(Article_length) + Article_poll_type + year + (1 | Article_ID), 
                  data=data_h4b, family=binomial(link="logit"), 
                  control=glmerControl(optimizer="bobyqa", 
                                       optCtrl=list(maxfun=2e5)))
summary(model504)

# Fixed effects specification
model504c <- glm.cluster(Change_correct ~ Broadsheet + TV + type_of_moe + + log(Article_length) + Article_poll_type + year, 
                         data=data_h4b, 
                         cluster='Article_ID',
                         family=binomial(link="logit"))
summary(model504c)




# Create texreg output tables ---------------------------------------------

# Models H1 - H3

screenreg(l=list(model100, model200, model300))
htmlreg(l=list(model100, model200, model300),
        custom.coef.names=c(NA, 
                            "Pollster WAPOR index", 
                            "Article length (log)", 
                            "Full article (Ref.: In-line reference)", 
                            "Year = 2012",
                            "Year = 2017",
                            "Margin of error in survey report", 
                            "Traditional MOE", 
                            "Uncertainty interval"),
        reorder.coef =c(1,2,7,8,9,3,4,5,6),
        groups = list("Type of MOE in survey report (Ref.: None)"=4:5, 
                      "Control variables"=6:9),
        custom.note = "%stars. The dependent variable concerns the quality of the news report. Model 1 is a linear regression with random effects for pollster. Models 2 and 3 are binary logistic models with random effects for pollster.",
        caption="",
        custom.model.names = c("(1) WAPOR Index", 
                               "(2) Margin of Error", 
                               "(3) Margin of Error"),
        file="models_h1_3.doc")



# Models H4

screenreg(l=list(model400, model451))
htmlreg(l=list(model400, model451),
        custom.coef.names=c(NA, 
                            "Traditional MOE", 
                            "Uncertainty interval",
                            "Article length (log)", 
                            "In-line poll reference (Ref.: full article)",
                            "Year = 2012",
                            "Year = 2017"),
        groups = list("Type of MOE in survey report (Ref.: None)"=2:3, 
                      "Control variables"=4:7),
        custom.note = "%stars. The dependent variable concerns the quality of the news report. Binary logistic regression models with random effects for article.",
        caption="",
        custom.model.names = c("(1) Difference Between Parties", 
                               "(2) Change Over Time"),
        file="models_h4.doc")


# Models H5

screenreg(l=list(model501, model502, model503, model504))
htmlreg(l=list(model501, model502, model503, model504),
        custom.coef.names=c(NA,
                            "Broadsheet newspaper",
                            "TV news website",
                            "Pollster WAPOR index", 
                            "Article length (log)", 
                            "Full article (Ref.: in-line reference)", 
                            "Year = 2012",
                            "Year = 2017",
                            "Traditional MOE", 
                            "Uncertainty interval"),
        reorder.coef =c(1,2,3,4,9,10,5,6,7,8),
        groups = list("Media type (Ref.: Tabloid newspaper)"=2:3,
                      "Type of MOE in survey report (Ref.: None)"=5:6, 
                      "Control variables"=7:10),
        custom.note = "%stars. The dependent variable concerns the quality of the news report. Model 1 is a linear regression with random effects for pollster. Model 2 is a binary logistic models with random effects for pollster. Models 3 and 4 are binary logistic models with random effects for article.",
        caption="",
        custom.model.names = c("(1) WAPOR Index", 
                               "(2) Margin of Error", 
                               "(3) Differences Between Parties",
                               "(4) Change Over Time"),
        file="models_h5.doc")





# Fixed effects versions with clustered standard errors

# H1- H3
screenreg(l=list(model100c, model200c, model300c))
htmlreg(l=list(model100c, model200c, model300c),
        custom.coef.names=c(NA, 
                            "Pollster WAPOR index", 
                            "Article length (log)", 
                            "Full article (Ref.: In-line reference)", 
                            "Year = 2012",
                            "Year = 2017",
                            "Margin of error in survey report", 
                            "Traditional MOE", 
                            "Uncertainty interval"),
        reorder.coef =c(1,2,7,8,9,3,4,5,6),
        groups = list("Type of MOE in survey report (Ref.: None)"=4:5, 
                      "Control variables"=6:9),
        custom.note = "%stars. The dependent variable concerns the quality of the news report. Model 1 is a linear regression with errors clustered by poll company. Models 2 and 3 are binary logistic models with errors clustered on poll company.",
        caption="",
        custom.model.names = c("(1) WAPOR Index", 
                               "(2) Margin of Error", 
                               "(3) Margin of Error"),
        file="models_fixef_h1_3.doc")

# H4
screenreg(l=list(model400c, model451c))
htmlreg(l=list(model400c, model451c),
        custom.coef.names=c(NA, 
                            "Traditional MOE", 
                            "Uncertainty interval",
                            "Article length (log)", 
                            "In-line poll reference (Ref.: full article)",
                            "Year = 2012",
                            "Year = 2017"),
        groups = list("Type of MOE in survey report (Ref.: None)"=2:3, 
                      "Control variables"=4:7),
        custom.note = "%stars. The dependent variable concerns the quality of the news report. Binary logistic regression models with errors clustered on poll company.",
        caption="",
        custom.model.names = c("(1) Difference Between Parties", 
                               "(2) Change Over Time"),
        file="models_fixef_h4.doc")


# H5

screenreg(l=list(model501c, model502c, model503c, model504c))
htmlreg(l=list(model501c, model502c, model503c, model504c),
        custom.coef.names=c(NA,
                            "Broadsheet newspaper",
                            "TV news website",
                            "Pollster WAPOR index", 
                            "Article length (log)", 
                            "Full article (Ref.: in-line reference)", 
                            "Year = 2012",
                            "Year = 2017",
                            "Traditional MOE", 
                            "Uncertainty interval"),
        reorder.coef =c(1,2,3,4,9,10,5,6,7,8),
        groups = list("Media type (Ref.: Tabloid newspaper)"=2:3,
                      "Type of MOE in survey report (Ref.: None)"=5:6, 
                      "Control variables"=7:10),
        custom.note = "%stars. The dependent variable concerns the quality of the news report. Model 1 is a linear regression with standard errors clustered by pollster. Model 2 is a binary logistic models with standard errors clustered by pollster. Models 3 and 4 are binary logistic models with standard errors clustered by article.",
        caption="",
        custom.model.names = c("(1) WAPOR Index", 
                               "(2) Margin of Error", 
                               "(3) Differences Between Parties",
                               "(4) Change Over Time"),
        file="models_fixef_h5.doc")



# Additional models for WAPOR items

screenreg(l=list(model101, model102))
htmlreg(l=list(model101, model102),
        custom.coef.names=c(NA, 
                            "Sample size in survey report", 
                            "Article length (log)", 
                            "Full article (Ref.: In-line reference)", 
                            "Year = 2012",
                            "Year = 2017",
                            "Data collection method in survey report"),
        reorder.coef =c(1,2,7,3,4,5,6),
        groups = list("Control variables"=4:7),
        custom.note = "%stars. The dependent variable concerns the quality of the news report. Models are binary logistic models with random effects for pollster.",
        caption="",
        custom.model.names = c("(1) Sample size", 
                               "(2) Data collection method"),
        file="models_extra_wapor.doc")

# Information about the R session
sessionInfo()

